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procedure  sum  (x,  n,  m,  result,  fail); 
value  n,  m;  integer  n,  m;  real  result; 

array  x;  label  fail; 

t  i  ,  1 

begin  comment  This  Algol  60  procedure  is  an  implementation  of  the 
floating-point  summation  technique  described  in  Malcolm  (1971).  This 
implementation  is  machine-independent  in  the  sense  that  it  will  work  on 
any  computer  having  a  floating-point  number  system  F  characterized  as 
follows:  Each  number  x€F  has  a  radix-f3  t-digit  fraction  where  t  >  1  . 
The  radix  0  can  be  any  positive  integer  greater  than  1  .  The  exponent 
e  is  assumed  to  lie  in  the  range 
b  <  e  <  B  , 

where  b  <  0  and  B  >  t  .  Each  nonzero  x€F  has  the  representation 
x  —  +  .d^  ...  dfc-  0®  , 

where  d^  ...  ,  are  integers  satisfying 
0  <  dA  <  P-1  ,  (i-1,  ...  ,t)  . 

The  number  0  is  contained  in  F  ,  but  no  assumption  is  made  about  its 
representation.  All  floating-point  operations  (e.g.,  addition  and  multi¬ 
plication)  are  assumed  to  result  in  either  0  or  a  normalized  floating¬ 
point  number  contained  in  F  .  The  machine  may  do  either  proper  rounding 
or  chopping  (truncation).  (Note  that  this  definition  of  F  excludes 
machines  using  extra-length  accumulators  for  intermediate  arithmetic. 
However,  this  algorithm,  is  seldom  needed  on  such  machines.) 

The  parameters  *}  and  t  of  F  are  automatically  computed  at 
execution  time  by  a  technique  described  in  Malcolm  (1972).  Since  the 
range  of  the  floating-point  exponent  cannot  be  determined  automatically, 
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the  input  parar.eter  m  is  used  for  allocating  the  set  of  accumulators 
used  by  the  algorithm. 

Provided  no  overflow  or  underflow  occurs,  and  none  of  the  x[i] 

are  larger  than  10m  ,  or  smaller  than  10 “ra  ,  in  magnitude,  and 

n  <  p^+1/l6  ,  where  l  =  [t/2\  ,  then 
n 

result  «  E  x[i] 
i=l 

is  returned  with  nearly  full-precision  accuracy.  The  bound  on  the 
relative  error  is  given  by  Theorem  2  in  Malcolm  (l^fl)  as 

r(t+l)/Uog  I6j1  p1_t  . 

P 

If  any  of  the  x[i]  are  larger  than  10m  or  smaller  than  10~m  ,  then 
the  error  exit  fail  is  taken.  ; 

Boolean  rnd;  integer  beta,  t,  t2,  nu,  L,  U; 
procedure  ENVRON  (beta,  t,  rnd); 

Boolean  rnd;  integer  beta,  t; 

begin  comment  This  procedure  is  an  Algol  60  translation  of  the  (first) 
Fortran  subroutine  ENVRON  given  in  Malcolm  (l$Ff2).  ; 

real  a,  b,  e; 

for  e  :=  2,  2xe  while  (a+l)-a=l  do  a  :=  e; 

for  e  :=  2,  2xe  while  a+b=a  do  b  :=  e; 

beta  :=  (a+b)-a;  rnd  :=  a+(beta-l)  >  a;  t  :=  0; 
for  a  :=  1,  betaXa  while  (a+l)-a=l  do  t  :=  t+1 

end  ENVRON; 

ENVRON  (beta,  t,  rnd);  t2  :=  t-f 2  ;  nu  :=  4n(l6)/£n(beta) ; 

U  :=  entier  (mX£n(lO)/(jin(beta)xnu))  +  1; 

L  :=  entier((-mx-£n(lO)/jtn(beta)  -  t2)/nu); 
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comment  In  the  notation  of  Malcolm  (lS^l),  t  =  t2  is  the  padding  that 


each  of  the  numbers  added  to  the  accumulators  will  have.  Each  of  the  x[i] 
will  be  split  into  two  helves  (i.e.  q=2)  having  the  last  t2  digits  equal 
zero.  The  variable  nu  above  is  used  for  v  defined  in  Equation  (2)  of 
Malcolm  (1971).  The  value  for  nu  computed  above  is  rather  arbitrary 
and  was  chosen  to  make  nu  sufficiently  smaller  than  t2  .  The  variables 
U  and  L  are  the  upper  and  lower  bounds  on  the  indices  of  the  accumulators 
which  are  declared  in  the  following  block.  They  are  chosen  to  allow  the 
x[i]  to  range  from  10"m  to  10m  in  magnitude.  In  slightly  different 
notation,  they  are 

U  =  rV(vxfog10p)l  , 

L  =  l(-y*°g10P  -  L  t/2J )/ vj  ; 

begin  array  accumulators[L:U] ;  integer  ex; 
real  xL,  xH; 


integer  procedure  e(x); 
value  x;  real  x; 

begin  comment  This  procedure  computes  the  exponent  e  of  the 
floating-point  number  x  .  ; 

real  y,  q;  integer  ex; 
x  :=  abs(x);  ex  :=  0;  for  y  :=  l,q 
while  x>y  do  begin  ex  :=  ex+1;  q  :=  betaxy;  end; 
for  y  :=  q,  y/beta  while  x<y  do  ex  :=  ex-1; 


e  :  = 
end  e; 
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comment  initialize  the  array  of  accumulators} 


for  i*=L  step  1  until  U  do  accumulators! i]  :=  0; 

comment  accumulate  the  nonzero  x[i]a; 
for  i  :=1  step  1  until  n  do  if  x[i]/0  then 

begin  ex  •=  e(x[i] ); 

if  entier(ex/nu)>U  V  ex-t2<Lxnu  then  go  to  fail; 
comment  Now  the  x[i]  is  split  into  a  high-  and  low-order 
part,  xH  and  xL.  The  method  used  here  is  to  add  the  proper 
power  of  P  to  x[i]  to  force  it  to  preshift  t2  digits 
to  the  right  and  then  either  truncate  or  round  the  last  t2 
significant  digits.  Then  the  same  power  of  P  is  subtracted 
to  cause  a  post  normalization  which  brings  in  t2  trailing 
zero  digits.  The  resulting  high-order  part  of  x[i]  is  then 
subtracted  from  x[i]  to  produce  the  low-order  part  such 
that  the  sum  of  the  high-  and  low-  order  parts  is  exactly 
equal  to  x[i].  This  method  of  splitting  a  floating-point 
number  into  two  halves  is  similar  to  that  given  by  Dekker 

(1971).  ; 

xH  :=  betat(ex-l+t2);  xH  :=  (xH+x[i])  -  xH; 
xL  :=  x[i]  -  xH; 

comment  xH  and  xL  can  now  be  added  to  the  appropriate 
accumulators.  ; 

accumulators[entier(ex/nu)]  :=  xH; 
accumulators[entier((ex-t2)/nu)]  :=  xL 
end;  comment  Now  sum  the  accumulators  in  decreasing  order.  ; 
result  :=  0; 

for  i  :=U  step  -1  until  L  do 
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result  ja  result  +  accumulators!  i] 


end 


end 


sum 
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